A partir da base de dados dos preços de combustíveis na cidade de Goiânia, será feita a filtragem listada acima. Foi selecionado um único posto revendedor - o posto com maior número de observações no período citado - contendo 365 valores divididos em 188 datas distintas; a seguir, foi tomada a média de cada mês essa será a variável estudada; foi feito a imputação de valores a cinco meses que não possuíam dados (nov-2018, dez-2018, fev-2019, set-2020, fev-2023) por meio da média dos seus respectivos anos; por fim,
# Retirando 'Revenda' e 'Produto'.
dados = dados3[, -c(1,2)]
summary(dados)
## Data.da.Coleta Valor.de.Venda mes ano
## Min. :2018-01-03 Min. :3.699 Length:365 Length:365
## 1st Qu.:2019-12-04 1st Qu.:4.690 Class :character Class :character
## Median :2021-04-07 Median :4.990 Mode :character Mode :character
## Mean :2021-03-13 Mean :5.396
## 3rd Qu.:2022-09-29 3rd Qu.:5.890
## Max. :2023-12-26 Max. :7.870
# Meses faltantes
m <- check_missing_months(dados, "Data.da.Coleta")
## [1] 5
print(m)
## [1] "2018-11" "2018-12" "2019-02" "2020-09" "2023-02"
# xodo: 5 faltas "2018-11" "2018-12" "2019-02" "2020-09" "2023-02"
# hiper moreira: 18 faltas
# viena: 5
# Monaco: 4
Como há 5 meses sem nenhuma observação foi utilizado da imputação de dados. Para cada mês listado acima será imputado o valor da média do respectivo ano, portanto, novembro e dezembro de 2018 terão o mesmo valor.
# Calculando médias anuais
aux = dados %>% group_by(ano) %>% mutate(media_anual = mean(Valor.de.Venda))
#unique(aux[, c(4, 5)])
# "2018-11" "2018-12" "2019-02" "2020-09" "2023-02"
# Adicionar os seguintes valores
linha1 = data.frame(matrix(c("2018-11-01", 4.661141, "2018-11", "2018",
"2018-12-01", 4.661141, "2018-12", "2018",
"2019-02-01", 4.665063, "2019-02", "2019",
"2020-09-01", 4.460714, "2020-09", "2020",
"2023-02-01", 5.702877, "2023-02", "2023"),
byrow= T, ncol=4))
colnames(linha1) <- colnames(dados)
linha1[, 1] <- as.Date.character(linha1[, 1])
linha1[, 2] <- as.numeric(linha1[, 2])
d1 = unique(dados_imput[, c(3, 5)])
a = qcc::qcc(d1$media_mensal, type = 'xbar.one',
title = 'Carta da média amostral (Dados finais)')
a = qcc::ewma(d1$media_mensal, title = 'Carta MMEP (Dados finais)')
a = qcc::cusum(d1$media_mensal, title = 'Carta "Soma cumulativa" (Dados finais)')
# AGRUPAMENTO TRIMESTRAL
m1 = matrix(d1$media_mensal, ncol = 3, byrow = T)
m11 = matrix(d1$mes, ncol=3, byrow = T)
a = qcc::qcc(data = m1, type = "S")
a = qcc::qcc(data = m1, type = "R")
A aplicação das cartas aos dados autocorrelacionados indica uma quantidade muito grande de pontos fora de controle. É de se esperar que o processo não seja bem medido e tais dados façam com que as cartas tornem-se, frequentemente, errôneas.
par(bg='gray')
d1 = unique(dados_imput[, c(3, 5)])
serie = ts(d1[, 2], start=c(2018, 1), frequency = 12)
#serie = ts(dados_imput$media_mensal, start=c(2018, 1), end=c(2023, 12), frequency = 12)
plot(serie, ylab="Valor", xlab="Obs.")
forecast::autoplot(decompose(serie))
Verificando a presença de autocorrelação entre as observações da série nota-se que há o decaimento exponencial à medida que se aumenta o lag.
par(bg='gray')
acf(d1$media_mensal, main = "Correlação entre as médias mensais")
Após toda organização do conjunto de dados, foram testados diversos modelos SARIMA. A seguir há uma relação dos modelos que apresentaram resíduos com: normalidade, estacionariedade e independência.
Outros modelos que também apresentaram todos pressupostos.
#aux = forecast::auto.arima(serie)
ar_ = forecast::Arima(serie, order = c(2, 3, 0),
seasonal = c(0,0,0), method = "ML")
summary(ar_)
## Series: serie
## ARIMA(2,3,0)
##
## Coefficients:
## ar1 ar2
## -0.9819 -0.4259
## s.e. 0.1076 0.1073
##
## sigma^2 = 0.2455: log likelihood = -48.96
## AIC=103.93 AICc=104.3 BIC=110.63
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.001505884 0.4779935 0.353291 0.1276067 6.67558 0.3425822
## ACF1
## Training set -0.09204694
# AMBOS PARÂMETROS DO MODELO ARIMA SÃO NEGATIVOS!!
shapiro.test(residuals(ar_))
##
## Shapiro-Wilk normality test
##
## data: residuals(ar_)
## W = 0.96813, p-value = 0.06461
tseries::adf.test(residuals(ar_))
## Warning in tseries::adf.test(residuals(ar_)): p-value smaller than printed
## p-value
##
## Augmented Dickey-Fuller Test
##
## data: residuals(ar_)
## Dickey-Fuller = -6.1038, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary
Box.test(residuals(ar_), type="Ljung")
##
## Box-Ljung test
##
## data: residuals(ar_)
## X-squared = 0.63581, df = 1, p-value = 0.4252
forecast::checkresiduals(ar_)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(2,3,0)
## Q* = 23.918, df = 12, p-value = 0.02087
##
## Model df: 2. Total lags used: 14
#agrup = qcc::qcc.groups(as.vector(ar_$residuals),
# as.vector(unique(dados_imput[, 3])))
g1 = qcc::qcc(ar_$residuals, type = "xbar.one", plot = T,
xlab="Grupo (Mês)",
ylab="Estatística do grupo",
add.stats = F,
label.limits = c("LCI", "LCS"))
g1$violations$beyond.limits
## [1] 57
g2 = qcc::ewma(ar_$residuals)
g3 = qcc::cusum(ar_$residuals)
#plot(g1, restore.par = F)
#head(d1)
d1[order(d1$media_mensal, decreasing = T), ]
## # A tibble: 72 × 2
## # Groups: mes [72]
## mes media_mensal
## <chr> <dbl>
## 1 2022-05 7.63
## 2 2022-04 7.62
## 3 2022-06 7.61
## 4 2022-03 7.53
## 5 2022-01 7.36
## 6 2021-11 7.34
## 7 2021-12 7.03
## 8 2021-10 7.02
## 9 2022-02 6.99
## 10 2021-09 6.66
## # ℹ 62 more rows
d1[56:58, ]
## # A tibble: 3 × 2
## # Groups: mes [3]
## mes media_mensal
## <chr> <dbl>
## 1 2022-08 5.10
## 2 2022-09 4.96
## 3 2022-10 4.99
g1$statistics[56:58]
## 56 57 58
## 0.9355005 1.6073610 0.3819759
g1 = qcc::qcc(ar_$residuals, type = "xbar.one", plot = F,
add.stats = F)
# 1. Convert qcc object to a data frame for ggplot2
qcc_data <- data.frame(
x = 1:length(g1$data), # x-axis (sample number)
Value = g1$data, # y-axis (data values)
CL = mean(g1$data[, 1]), # Center Line
UCL = g1$limits[2], # Upper Control Limit
LCL = g1$limits[1], # Lower Control Limit
violations = ifelse(g1$data > g1$limits[2] | g1$data < (g1$limits[1] - 3*(g1$limits[2] - g1$limits[1])), "Violação", "Sob controle"))
# 2. Create the ggplot2 plot
gg1 <- ggplot(qcc_data, aes(x = x, y = Value)) +
theme_classic() +
geom_point(aes(color = violations)) + # Points, colored by violations
geom_line() + # Connect the points
geom_hline(yintercept = qcc_data$CL[1], linetype = "solid", color = "black",
linewidth = .1) + # Center line
geom_hline(yintercept = qcc_data$UCL[1], linetype = "dashed", color = "red",
linewidth = .3) + # Upper control limit
geom_hline(yintercept = qcc_data$LCL[1], linetype = "dashed", color = "blue",
linewidth = .3) + # Lower control limit
labs(title = "", x = "Grupo (Mês)", y = "Estatística do grupo") + # Labels
scale_color_manual(values = c("Violação" = "red", "Sob controle" = "black")) + # Custom colors
annotate(geom = 'text', label = c(paste0("LSC = ", round(qcc_data$UCL[1], 3)),
paste0("LC = ", round(qcc_data$CL[1], 3)),
paste0("LIC = ", round(qcc_data$LCL[1], 3))),
size = 3, x = c(5,1,5), y = c(qcc_data$UCL[1] + .1 ,
qcc_data$CL[1] + .1,
qcc_data$LCL[1] + .1)) +
theme(legend.position = "none") # Remove the legend
plotly::ggplotly(gg1)
Gráficos utilizados no relatório.
ar1 = matrix(ar_$residuals, ncol = 3, byrow = T)
ar11 = matrix(d1$mes, ncol = 3, byrow = T)
ar111 = matrix(c(1:72), ncol = 3, byrow = T)
am1 = qcc::qcc(ar1, type = "R")
g1$violations$beyond.limits
## [1] 57
#ar11[g1$violations$beyond.limits, ]
#qcc::qcc(ar1, type="S")
ar1 = matrix(ar_$residuals, ncol = 4, byrow = T)
ar11 = matrix(d1$mes, ncol = 4, byrow = T)
ar111 = matrix(c(1:72), ncol = 4, byrow = T)
am1 = qcc::qcc(ar1, type = "R")
am1$violations$beyond.limits
## [1] 14 15
#ar11[am1$violations$beyond.limits, ]
qcc::qcc(ar1, type="S")
## List of 11
## $ call : language qcc::qcc(data = ar1, type = "S")
## $ type : chr "S"
## $ data.name : chr "ar1"
## $ data : num [1:18, 1:4] 0.00104 0.10152 0.28529 0.16085 0.04526 ...
## ..- attr(*, "dimnames")=List of 2
## $ statistics: Named num [1:18] 0.0244 0.4928 0.363 0.2921 0.2695 ...
## ..- attr(*, "names")= chr [1:18] "1" "2" "3" "4" ...
## $ sizes : int [1:18] 4 4 4 4 4 4 4 4 4 4 ...
## $ center : num 0.457
## $ std.dev : num 0.496
## $ nsigmas : num 3
## $ limits : num [1, 1:2] 0 1.04
## ..- attr(*, "dimnames")=List of 2
## $ violations:List of 2
## - attr(*, "class")= chr "qcc"
ar1 = matrix(ar_$residuals, ncol = 3, byrow = T)
ar11 = matrix(d1$mes, ncol = 3, byrow = T)
ar111 = matrix(c(1:72), ncol = 3, byrow = T)
am1 = qcc::qcc(ar1, type = "S")
g1$violations$beyond.limits
## [1] 57
#ar11[g1$violations$beyond.limits, ]
# 1. Convert qcc object to a data frame for ggplot2
qcc_data <- data.frame(
x = 1:length(am1$statistics), # x-axis (sample number)
Value = am1$statistics, # y-axis (data values)
CL = mean(am1$statistics), # Center Line
UCL = am1$limits[2], # Upper Control Limit
LCL = am1$limits[1], # Lower Control Limit
violations = ifelse(am1$statistics > am1$limits[2] | am1$statistics < (am1$limits[1] - 3*(am1$limits[2] - am1$limits[1])), "Violação", "Sob controle"))
# 2. Create the ggplot2 plot
gg1 <- ggplot(qcc_data, aes(x = x, y = Value)) +
theme_classic() +
geom_point(aes(color = violations)) + # Points, colored by violations
geom_line() + # Connect the points
geom_hline(yintercept = qcc_data$CL[1], linetype = "solid", color = "black",
linewidth = .1) + # Center line
geom_hline(yintercept = qcc_data$UCL[1], linetype = "dashed", color = "red",
linewidth = .3) + # Upper control limit
geom_hline(yintercept = qcc_data$LCL[1], linetype = "dashed", color = "blue",
linewidth = .3) + # Lower control limit
labs(title = "", x = "Grupo (Mês)", y = "Estatística do grupo") + # Labels
scale_color_manual(values = c("Violação" = "red", "Sob controle" = "black")) + # Custom colors
annotate(geom = 'text', label = c(paste0("LSC = ", round(qcc_data$UCL[1], 3)),
paste0("LC = ", round(qcc_data$CL[1], 3)),
paste0("LIC = ", round(qcc_data$LCL[1], 3))),
size = 3, x = c(5,1,5), y = c(qcc_data$UCL[1] + .1 ,
qcc_data$CL[1] + .1,
qcc_data$LCL[1] + .1)) +
theme(legend.position = "none") # Remove the legend
plotly::ggplotly(gg1)